1. Field of the Invention
This invention relates generally to a system for analyzing signals and data and, more specifically but without limitation, to a system for analyzing non-linear and/or non-stationary signals and data.
2. Description of Related Art
Data analysis is an essential part of pure and applied research in many disciplines. In many practical applications, the raw data have characteristics that change over time, commonly referred to as being non-stationary and, additionally, very often result from non-linear processes. Unfortunately, the majority of existing methods for analyzing data are designed to treat stationary, linear data.
Another common and serious problem of data analysis is the existence of noise and/or non-stationary trend. Common practice to deal with these problems has involved application of band pass filters to the data. However, these filters are Fourier based and, as such, typically result in the introduction of spurious harmonics in non-stationary data. Therefore, Fourier-based filters have limited utility and practical value for use with non-stationary and/or non-linear data. In addition, the (low frequency) signal trend can carry significant and useful information about the process being analyzed, and thus should not be simply filtered out. Prior art methods for processing non-stationary data include Fourier analysis, wavelet analysis, the Wigner-Ville distribution, the evolutionary spectrum, the empirical orthogonal function expansion, and the empirical mode decomposition. These prior art methods can be briefly described as follows.
Fourier Analysis.
Traditional methods of time-frequency-energy analysis are based on Fourier transforms and are designed for stationary and linear processes. The application of these methods to analysis of non-stationary, non-linear data can give misleading results. For example, the Fourier spectrum defines uniform harmonic components globally. Therefore, Fourier analysis needs many additional harmonic components in order to simulate non-stationary data, which are not uniform globally. As a result, Fourier analysis tends to spread signal energy over a wide frequency range. As is well known by those having skill in the art, the faster the change in the time domain, the wider the frequency range. Unfortunately, many of the components, that are added in order to simulate the non-stationary nature of the data in the time domain, divert energy to a much wider frequency domain. For example, a single impulse, a signal whose deviation from constancy occurs at a single moment in time, requires infinitely many frequencies with identical power to represent it. Constrained by energy conservation, the spurious harmonics that are added and the wide frequency spectrum required to simulate the non-linearity cannot faithfully represent the true energy density in the resulting frequency space.
Further, Fourier spectral analysis utilizes linear superposition of trigonometric functions. Such analysis needs additional harmonic components in order to simulate the effects of non-linearity, such as deformed wave profiles. Whenever the form of the data deviates from a pure sine or cosine function, the Fourier spectrum will contain harmonics. As explained above, both non-stationarity and non-linearity can induce spurious harmonic components that cause energy spreading. The resulting consequence is a misleading or incorrect time-frequency distribution of signal energy for non-linear and/or non-stationary data.
Many data analysis methods have been developed based on Fourier transforms. The spectrogram is the most basic and common method, which is a finite-time window Fourier spectral analysis that is repeated in moving-time windows. By successively sliding the window along the time axis, a time-frequency-energy distribution is obtained. Since such a distribution relies on the traditional Fourier spectral analysis, the method assumes the data to be piecewise stationary. This assumption is not valid for most non-stationary data. Even if the data were piecewise stationary, it is highly unlikely in most cases that the window size adopted would coincide with the stationary time scale. Furthermore, there are practical difficulties in applying the method. In order to localize an event in time with good temporal precision, the window width must be narrow. Unfortunately, the frequency resolution worsens as window size decreases. Although the conflict in these requirements can be mitigated by different techniques, it renders the applicability of Fourier analysis to non-linear, non-stationary data of limited use.
Wavelet Analysis.
Wavelet analysis, which has become extremely popular during the past decade, is an attempt to overcome the problems of windowed Fourier analysis by utilization of a basis of functions to represent a signal that contains elements having different time scales. This approach allows wavelet analysis to detect changes that occur rapidly, i.e., those on a small time scale, with good temporal resolution but poor frequency resolution, or slowly, i.e., those on a large time scale, with good frequency resolution but poor temporal resolution. More specifically, the wavelet analysis approach is essentially an adjustable-window Fourier spectral analysis. Wavelet analysis is useful for analyzing data with gradual frequency changes. Primary applications of wavelet analysis have been in areas of edge detection and audio and image signal compression. Limited applications also include analysis of time-frequency distribution of energy in time series and of two-dimensional images such as fingerprints. Unfortunately, the finite length of the basic wavelet function results in energy leakage across different levels of resolution in a multi-resolution analysis, which causes quantitative applications of the time-frequency-energy distribution to be more difficult.
Sometimes, the interpretation of the wavelet can also be counterintuitive. For example, the more temporally localized the basic wavelet, the higher the frequency range will be. Therefore, to define a change occurring locally, the analytic result may occur in a high frequency range. In other words, if a local event occurs only in a low frequency range, the effects of that local event may only appear in a high frequency range. In many applications, the interpretation of such a result would be difficult if not impossible.
Another difficulty with wavelet analysis is that it is not adaptive. Once the basic wavelet is selected, the basis for the analysis is completely determined to the extent obtainable from the selected basic wavelet, and all information of the input signals is represented in terms of that basis. Although the basis can be specially selected for an individual application, the information obtained depends heavily on the properties inherent to that basis rather than solely on the intrinsic properties of the signals being studied. Malvar wavelets, Wavelet Packets, and Matching Pursuit methods have been developed to overcome some of these limitations to more accurately represent a signal having dynamics that vary with time and that include both stationary and non-stationary characteristics. Unfortunately, these developments in wavelet analysis continue to suffer from the representation of the signal information in terms of a pre-selected basis of functions that often has little or nothing to do with the dynamics and other characteristics of the input signals.
The Wigner-Ville Distribution.
The Wigner-Ville distribution, sometimes referred to as the Heisenberg wavelet, is the Fourier transform of the central covariance function. The difficulty with this method is the severe cross terms indicated by the existence of negative power for some frequency ranges. The Wigner-Ville distribution has been used to define wave packets that reduce a complicated data set into a finite number of simple components. Although this approach can be applied to a wide variety of problems, applications to nonstationary or nonlinear data are extremely complicated. Further, such applications again suffer from the same limitation of the other prior art methods described above in that the bases for representation of information are not derived from the data itself.
Evolutionary Spectrum.
The Evolutionary Spectrum approach is used to extend the classic Fourier spectral analysis to a more generalized basis, namely from sine or cosine functions to a family of orthogonal functions indexed by time and defined for all real frequencies. Thus, the original signal can be expanded into a family of amplitude modulated trigonometric functions. The problem with this approach, which severely limits its applicability, is due to the lack of means for adequately defining the basis. In principle, the basis has to be defined a posteriori in order for this method to work. To date, no systematic method of constructing such a basis exists. Therefore, it is impossible to construct an evolutionary spectrum from a given set of data. As a result, applications of the evolutionary spectrum method have changed the approach from data analysis to data simulation. Thus, application of the evolutionary spectrum approach involves assumptions causing the input signal to be reconstituted based on an assumed spectrum. Although there may be some general resemblance of a simulated input signal to the corresponding real data, it is not the data that generated the spectrum. Consequently, evolutionary spectrum analysis has very little useful application, if any.
The Empirical Orthogonal Function Expansion.
Empirical orthogonal function expansion (“EOF”), also known as “principal component analysis” or “singular value decomposition,” provides a means for representing any function of state and time as a weighted sum of empirical eigenfunctions that form an orthonormal basis. The weights are allowed to vary with time. EOF differs from the other methods described hereinabove in that the expansion basis is derived from the data. The critical flaw of EOF is that it only gives a variance distribution of the modes defined by the basis functions, and this distribution by itself gives no indication of time scales or frequency content of the signal. In addition, any single component of the non-unique decomposition, even if the basis is orthogonal, does not usually provide any physical meaning. The Singular Spectral Analysis (“SSA”) method, which is a variation of EOF, is simply the Fourier transform of EOF. Unfortunately, since EOF components from a non-linear and non-stationary data set are not linear and stationary, use of Fourier spectral analysis generated by SSA is flawed. Consequently, SSA is not a genuine improvement of EOF. Although adaptive in nature, the EOF and SSA methods have limited applicability.
The Empirical Mode Decomposition.
The empirical mode decomposition (“EMD”) method, involves two major steps. The first step is the application of an algorithm for decomposing physical signals, including those that may be non-linear or non-stationary, into a collection of Intrinsic Mode Functions (“IMFs”), which are supposedly indicative of intrinsic oscillatory modes in the physical signals. More specifically, the cornerstone of the EMD method is the extraction of a signal baseline from a physical signal wherein the baseline is computed as the mean value of the upper and lower envelopes of the physical signal. The upper envelope is defined by cubic splines connecting the local maxima of the physical signal, and the lower envelope is defined by cubic splines connecting the local minima of the physical signal. The signal baseline is then extracted, or subtracted, from the original signal to obtain an IMF having the first and highest frequency present in the signal.
A goal of the first step is to obtain a well-behaved IMF prior to performing the second step: applying a Hilbert Transform to the IMF. “Well-behaved” means that the IMF should be a “proper rotation,” i.e., all local maxima are strictly positive and all local minima are strictly negative. This does not necessarily happen with one step of EMD, and thus a laborious, iterative “sifting” process is applied to the signal baseline and is terminated when a set of “stopping criteria” are satisfied, such as when the resulting IMF either becomes a proper rotation or when some other criteria are reached (e.g., a predetermined number of iterations exhausted) without obtaining a proper rotation. The stopping criterion is based on a combination of (i) limiting the amount of computational energy expended, and (ii) having the constructed IMF closely approximate the desired property. When the first IMF function has been obtained, it is subtracted from the signal and the process is repeated on the resulting lower frequency signal. This process is repeated again and again until the decomposition is completed in the window of signal being analyzed.
The sifting process has two goals: (i) to separate out high frequency, small amplitude waves that are “riding” atop, or superimposed on, larger amplitude, lower frequency waves, and (ii) to smooth out uneven amplitudes in the IMF being extracted. Unfortunately, these goals are often conflicting for non-stationary signals wherein riding waves may be isolated and/or are highly variable in amplitude. As a result, the sifting process must be applied cautiously as it can potentially obliterate the physically meaningful amplitude fluctuations of the original signal.
As noted, once the IMFs have been obtained, the second step of the EMD method is to apply the Hilbert Transform to the IMFs which, provided the IMFs are well-behaved, results in quantified instantaneous frequency measurements for each component as a function of time.
However, the EMD method suffers from a number of shortcomings that lead to inaccuracies in depiction of signal dynamics and misleading results in subsequent signal analysis, as follows:                (a) The construction of the IMF baseline as the mean of the cubic spline envelopes of the signal suffers from several limitations, including the time scale being defined only by the local extrema of the original signal, ignoring the locations of the rest of the critical points, such as inflection points and zero crossings, which are not preserved by the sifting process.        (b) The EMD transformation is window-based; the sifting procedure and other processing requires an entire window of data to be repeatedly processed. The resulting information for the entire window is not available until the window's processing has been completed. On average, this causes a delay in obtaining the resulting information of at least half the window length plus the average computational time for a window of data. As a result, an insurmountable problem is created for the EMD in real-time/online applications. In order to obtain the information quickly, the window length should be short. However, short windows have less accurate results due to boundary or edge effects, and are incapable of resolving frequency information on a time scale longer than the window length itself.        (c) The EMD method is computationally expensive, and also subject to uncertain computational expenses. The EMD method requires repeated sifting of components in order to obtain well-behaved IMFs or until stopping criteria are satisfied. Such a procedure may require numerous iterations and may not occur in finite time. Thus, the method often will not result in IMFs with the desired proper rotation property, even when sifting numerous times. Also, the IMFs are generally dependent on the parameters of the algorithm that define the stopping criteria for sifting.        (d) Overshoots and undershoots of the interpolating cubic splines generate spurious extrema, and shift or exaggerate the existing ones. Therefore, the envelope-mean method of the extraction of the baseline does not work well in certain cases, such as when there are gaps in the data or data are unevenly sampled, for example. Although this problem can be mitigated by using more sophisticated spline methods, such as the taut spline, such improvements are marginal. Moreover, splines are not necessarily well suited to approximate long timescale trends in real data.        (e) The EMD method does not accurately preserve precise temporal information in the input signal. The critical points of the IMFs, such as the extrema, inflection points, etc., are not the same as those of the original signal. Also, the EMD method, being exclusively determined by extrema, is deficient in its ability to extract weak signals embedded in stronger signals.        (f) Since the cubic splines can have wide swings at the ends, the envelope-mean method of the EMD method is particularly unsuitable for real-time applications or for applications utilizing a narrow window. The end effects also propagate to the interior and significantly corrupt the data, as the construction of IMFs progresses, as can be seen in FIGS. 1a-2c. FIG. 1a illustrates the application of the EMD to a test signal (top-most signal) to produce IMF components (displayed below test signal). This panel illustrates end effects, spline-related instabilities (most noticeable in bottom components), and inability to extract the readily apparent signal baseline. The intrinsic timescale decomposition, sometimes referred to herein as ITD, (shown in FIG. 1b) separates the same test signal (top-most signal) into stable components (displayed below test signal), demonstrating the fact that that ITD has no end effect propagation beyond the first two extrema at each level and allows correct identification of the trend (dashed line). FIG. 2a shows a brain wave input signal (electrocorticogram; ‘ECoG’) containing an epileptic seizure used to illustrate decomposition differences between EMD and ITD. EMD-generated (lower left panel) and ITD-generated (lower right panel) decompositions of the cumulative sum of the raw signal show that ITD, unlike EMD, does not generate extraneous components and correctly reveals large timescale variations of the signal (DC trend).        (g) The application of the Hilbert Transform to track instantaneous frequency is only appropriate when frequency is varying slowly with respect to amplitude fluctuations which, unfortunately is a condition not satisfied by many non-stationary signals. Moreover, proper rotations, which are not guaranteed by the EMD, are necessary for the existence of a meaningful, instantaneous frequency that is restrictive and local.        (h) Primarily due to the sifting procedure, the EMD method causes (i) a smearing of signal energy information across different decomposition levels, and (ii) an intra-level smoothing of energy and frequency information that may not reflect the characteristics of the signal being decomposed or analyzed. The potential negative effects of over-sifting include the obliteration of physically meaningful amplitude fluctuations in the original signal. Also, the inter-level smearing of energy and limitation of decomposition in a window of data create inaccuracies in the EMD's representation of the underlying signal trend, especially if the signal trend is longer than a single window.        (i) The EMD method cannot guarantee that the IMF components will be “proper rotations,” even with the sifting reiterations. Often, in practice, they are not.        (j) The EMD method is not well behaved if the data is unevenly sampled or if it is discontinuous, and, therefore, may not preserve intact phase and morphology characteristics of the signal.        
What is needed is a system for reliably analyzing non-linear and/or non-stationary signals and data capable of decomposing them and/or extracting precise time-frequency-energy distribution information.